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Abstract. In this paper, we present a parallel scheme to solve the pop- 
ulation balance equations based on the method of characteristics and the 
finite element discretization. The application of the method of character- 
istics transform the higher dimensional population balance equation into 
a series of lower dimensional convection-diffusion-reaction equations which 
can be solved in a parallel way. Some numerical results are presented to 
show the accuracy and efficiency. 



1. Introduction 

In this paper, we propose a parallel scheme to solve the population bal- 
ance equation (PBE) based on the application of the method of characteristics 
and the finite element method. The PBEs aries from the model of the in- 
dustrial crystallization process (see, e.g., HU [T2] and the reference cited 
therein). Recently, more and more researchers are interested in the numerical 
methods for PBEs (c.f. [U El H3 In PBEs, besides the normal space and 
time variables, the distribution of entities also depends on their own proper- 
ties which are referred to as internal coordinates. It is a high dimensional 
system of equations which is a big challenge from the computational point of 
view. In order to overcome this difficulty, we use the method of characteristics 
(c.f. [21 H]) to transfer the original problem to a series of lower-dimensional 
convection-diffusion-reaction problems which are defined on the characteris- 
tics curves and the spatial directions. Based on the data structure for the 
method of characteristics, a parallel implementation can be applied to do the 
simulation process that can improve the computational efficiency. 
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So far, there exists the alternating direction (operator splitting) method for 
the PBE by decomposing the original problem into two unsteady subproblems 
of smaller complexity (see, e.g., [UEJE])- In the two subproblems, the ordering 
of the data for the solution needs to be different since they are discretized in 
different direction (c.f. [1]). It is not so suitable for the parallel implementation 
and prevents the further improvement of the computation efficiency for the 



In the present paper, we use the method of characteristics to transform the 
PBE into a series of convection-diffusion-reaction equations on the character- 
istic curves in each time step. Then the finite element method is applied to 
solve the series of convection-diffusion-reaction problems. Furthermore, based 
on the data structure of the numerical scheme, a parallel scheme is constructed 
to solve the PBE based on the distributed memory. Some numerical results 
are provided to check the efficiency of this parallel method. 

The following of the paper will go as follows: Section [2] introduces the model 
problem under consideration and defines some notation. In Section [31 we de- 
scribe method of characteristics for solving the PBE. The finite element dis- 
cretization for the PBE is described in Section HI Then Section gives the 
parallel implementation way for the full discrete form of the PBE. The numer- 
ical results are given in Section [H] to validate the efficiency of the numerical 
method proposed in this paper. Some concluding remarks are given in the last 
section. 



Let f2 x be a simply connected domain in TZ d (d = 2 or 3) with Lipschitz 
continuous boundary <9f2 x , Qi = [£ m i n ,£ max ] C 1Z and T > 0. The state of the 
individual particle in the PBE equation may consists of the external coordi- 
nate x (x = (xi, ■ • • , Xd)), denoting its position in the physical space, and the 
internal coordinate £, representing the properties of particles, such as size, vol- 
ume, temperature etc.. A PBE for a solid process such as crystallization with 
one internal coordinate can be described by the following partial differential 
equation: 

Find z : (0, T] x Q e x fi x K such that 



| + G{£) | - eA^z + b(x) ■ V x z = f(t, £, x) in (0, T] x Q e x O x , 



where the diffusion coefficient e > is a given constant, A x and V x denote the 
Laplacian and gradient with respect to x, respectively, b is a given velocity and 
satisfies V x ■ b = 0, and / is a source function. Here G(t) > represents the 
growth rate of the particles that depends on £ but is independent of x and t. 



PBE. 



2. Model problem 



(2.1) 
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Furthermore, let us assume the data G(£), b, /, z- 1T ax an d ^bdry are sufficiently 
smooth functions for our error estimate analysis. 

Now we introduce some notation of the function spaces (c.f. [21 Let 
if m (f2 x ) denote the standard Sobolev space of functions with derivatives up 
to m in L 2 (Q X ) and the norm is defined by 



, n — \dx a 

0<|a|<m 



dx. 



where a denote a non- negative multi- index a = {a\, 
and 



a d}, \a\ = 



dxl 



«1 



• X 



ad 



We use (•, -) x and || ■ Hl^o*) to denote the L 2 -inner product and the associated 
norm in fl x , respectively, which are defined as follows 



[V,W] 



vwdx. and \\v\\ L 2 



(^) 



Let X be a Banach space with the norm || • Then we define 

C(Q,f,X) = (v : Cli — > X : v is continuous j, 

cPv 



W m '°°(ttf,X) = [vMi-tX: 
iy m '°°((0,T];X) = (v:(0,7]->X 



x 



dP 



< oo, < j < mj, 
< oo, < j < m j, 



where the derivatives cPv/dft and d^v/dP are understood in the sense of dis- 
tribution on D,£ and (0, T], respectively The norms in the above defined spaces 
are given as follows 



\\v\\c(n e] x) 

\\V \\W m '°°(n t ;X) 
\v\\w m <°°({a,T\;X) 



sup \\v 



max sup 

0<j<m eeSi 



X ■ 



x 



max sup 

0<J<m te ( ,T] 



OP 



X 



For spaces X, Y and Z, we use the short notation Z(Y(X)) := Z((0, T]; (Y(Qf, X)) 
in this paper. 



3. Method of characteristics 



In this section, we describe the method of characteristics (c.f. [21 lU E]) for 
the PBE ( 12. II) . The reason we adopt this method for the discretization in 
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the product space (0,T] x Q,? is that it has the suitable data structure for the 
parallel implementation which will be discussed in the following sections. 
First we set 

VKM) = (i + G(£) 2 ) 1/2 - 

Let the characteristic direction associated with the hyperbolic part of (12. ip . 
dz/dt + G(£)dz/d£, be denoted by s(t). Then 

(3 1) - 1 9 | G W d 

1 ' ' ds ipdt V df 

Then (12.11) can be written as 



V'ff - sA x z + b(x) • V x 2 = / in (0, T] x O f x tt y 



(3.2) 



z(0, £, x) = z init (£, x) in Sl e x Q y 

£ min , x) = 2 b dry(t, x) on (0, T] x fi x , 

z(t, £, x) = on (0, T] x ^ x dfi, 



We use uniform partitions for the time interval (0,T] and the internal coordi- 
nate interval fi^, respectively. Let r = T/iV, i = (£ max — £ min )/M, t n = nr, 
n = 0, 1, ■ • • ,N and £ m = £ min + mi, m = 0, 1, • • • , M. In order to satisfy the 
stability condition, we set 

(3-3) r < — - 

max^ min < £ < £max G(f) 

Then starting with z(0,£, x) = z init , £ m i n , x) = Zbd ry (£, x), the equation 
(13. 2p can be discreted in each sub-intervals (t n_1 ,t n ] x (d_i,£ m ] x fi x (n = 
1, 2, • • • , iV and m = 1, 2, • • • , M) as follows. 
First we compute 

(3.4) £ m = £ m -rG(£ m ). 

Actually, this is a first order discretization to obtain the approximation at 
the time level t = t n_1 for the following characteristic ordinary differential 
equation (c.f. [1]): 



(3-5) 



dt 

£{t n ) 



M - G{£) in [t n -\t n ), 



From the condition (13.31) . we have £ m >= £ m i n for m > 1. Then we compute 
the direction differential t/>|j at the node (t n ,£ m ) in the following way 



(3-6) 



^(t i £-mi x) ^(t , £mi x) 

T 

in— I „\ „,n „(+n—\ l) „\ i /i „,n x^l+n— 1 



where z(rX,x) := 0(t n_1 , C-i, x) + (1 - a n m )z{t n ~\ £ m , x) with < 
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In order to give the semi-discrete form of the PBE, we set z^(x) m z(t n ,£ m , x) 
Then the semi-discrete form of the PBE can be defined as follows: 



£ A x C(x) + b(x)V x C(x) = /£(x) in n x , 



(3.7) 



^m( x ) = ^it( x ) for a; e O x , 

^(x) = z bdry (r,x) for(0,T]xfi 
(x) = for m = 1, 2, ■ ■ • , M onflU 



where /»(x) = x), S£(x) = <CCi(x) + (1 - <CK -1 (*)- 

From the Taylor expansion method, we can derive the following error esti- 
mate for the semi-discrete form (\3.7\i 

(3.8) \\z(t n ,£ m ,x) -2;^(x)||c(c(x)) < Cr||^(t,^x)|| W ra,ao (w - 1 ,oo( X )), 
where the space X can be L 2 (fi x ) or if 1 (f2 x ). 

4. Finite element method 

In this section, we give the fully discrete form of the PBE by the finite 
element method. Let Vh be a finite element subspace of Hq(Q x ) which has the 
fc-th order accuracy (c.f. [2j [3] ) : 

(4.1) inf \\u - T; fc ||ffi (ni0 < Ch k \\u\\ Hm+ i inx) Vu G H m+ \Q X ). 



and 



(4.2) inf \\u - v h \\ L2{Ux) <Ch k+1 \\u\\ H m +Hnx) VueH m+1 (n x ). 

Based on the finite element space Vh, we can define the fully discrete form 
for the PBE as follows: 

For the n-th time step t = t n and m = 0, 1, • • ■ , M, find z^ nh E Vh such that 
(4-3) 

( m ' h T m * ,v h ) + a(z^ h ,v h ) = (&(x),v h ) Vv h G V h , 

a Q (zl l h , v h ) = a (z init (£ m , x),v h ) Wv h e V h , m = 1, • • • , M, 

a (z% h , v h ) = a (z hdry (t n , x),v h ) Wv h G V h , 

where z^ nh = a^^-i h + 0- ~ a m) z rn~h w ^ n a m being defined in Section [3] and 



/ (eVu ■ Vv + b(x) ■ Vit v)dx, 



a{u, v) 

a (u,v) = I Vu ■ Vucix. 

From the standard error estimate theory of the finite element method (c.f. 
[2j [3] ) , the fully discrete form (14. 3 j) has the following error estimates 

(4.4>nax : \\z(T,£ m ,x) - z% h \\ H i {Ux) < C(r + /i fc )||«|| w a,oo (wl ,ao ( ^*+i (ax))) 

Km<M 
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and 

(4.S^ax^||^(T,^ m ,x) -z^ h \\ L 2 {nx) < C(t + h k+1 )\\z\\ W 2, x{w i, x{H k+i {nx))) . 

5. A PARALLEL WAY 

In this section, we present a parallel scheme to solve the PBE ( 12. ip based on 
the full discrete (14.31) . Fortunately, from (14. 3p . we can find the finite element 
equation is independent for each m in any time step t n . Based on this property, 
we can construct a type of parallel scheme to implement the full discretization 
of the fully discrete PBE (Oil . 

Assume we use P processors to compute the PBE. Decompose the set 
{0, 1, 2, • ■ ■ , M} into P subsets mi, m 2 , ■ ■ • , rrip such that mi = {m = 0, 1, ■ ■ • , mi — 
1}, m p = {m p _i, m p _i + 1, - • • , m p — 1} (p = 2, ■ • • , P - 1) and m P = 
{mp-i, ■ ■ ■ ,mp — 1 = M}. In the p-th processor, the equation f!4.3[) is 
solved on the sub-intervals (t n ~ l ,t n ] x (^pCp-i] x Q x (n — 1,2, ■ ■ ■ , N, 
p = 1, 2, • • • ,P, £o = imin and £m = 4nax)- Because the growth rate of the 
particles G(£) > 0, the dependence of each point £ m is on the left (£ < £ m ), the 
solution z^ ~_i h in ^ ne P"th processor as the initial condition for the p + 1-th 
processor computing at the t n time step. 

We allocate the memory in the p-th processor (p = 1, • • • , P) to save the 
solutions -i,h> ' ' ' ' z m p -i,h an d the p-th processor {p = 1, • • • , P — 1) should 
send its saved solutions to the next p + 1-th processor after each time step 
computation. Obviously, for p — 1, we need to use the boundary condition 
-2bdry(A x ) and for p = P, the sending of solutions is not required since it is the 
last processor. Based on this distribution of the memory and the computation 
of the scheme (14. 3p . we can construct the following parallel algorithm for the 
PBE. 

Algorithm 5.1. Parallel algorithm for PBE 
For n = 1, 2, • • • , JV Do 

(1) On each processor, compute the solution z^ nh for m 6 m p (p = 1, 2, ■ ■ • , P) 

in this sub-interval (t n ~ 1 ,t n ] x {i mp -i , 4n p -i] ■ 

(2) For p = 1, 2, • • • , P— 1, send the solutions in the p-th processor z^ nh {m e 
m p ) to the p + 1-th processor. 

(3) If n < N, set n := n + 1 and go to Step 1. Else stop. 

6. Numerical results 

In this section, we provide some numerical results to validate the numerical 
scheme proposed in this paper. Let fi x = [0, 1] x [0, 1], = [0, 1], T = 1, e — 1 
and b(x) = (1, 1) T . We chose the functions f(t,£, x), z init (£, x) and Zbdry(£, x) 
such that the exact solution is 

z(t, £, x, y) = e~ at sin(7r£) sin(7rx) sin(7n/) 

with a = 0.1. The growth rate of the particles is G(£) = | + 2(1 — £)£. 
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First, we check the convergence order for the error estimates 
||e|| = max ||z(T,£,x) -2™ J L2(nx) 

l<m<M 

and 

Hi = max \\z(T,£,x) - z^ h \\ H x {nx) . 

l<m<M 

The convergence order of the linear and quadratic finite element method for 
the discretization in f2 x is shown in Tables [1] and [2j From Tables [1] and [2], we 
can find the numerical method with the linear and quadratic finite element 
method in the space direction has the reasonable convergence order. 

Table 1. Error and rate of convergence in the space direction 
for PI with t = l = h 2 



mesh 


size h 


IHIo 




IMIi 








error 


order 


error 


order 


2 


-2 


4.5702E-01 




2.6897E-00 




2 


-3 


1.4872E-01 


1.6197 


1.5128E-00 


0.8302 


2 


-4 


4.0481E-02 


1.8773 


7.8083E-01 


0.9541 


2 


-5 


1.0318E-02 


1.9721 


3.9359E-01 


0.9883 


2 


-6 


2.7230E-03 


1.9219 


1.9720E-01 


0.9970 



Table 2. Error and rate of convergence in the space direction 
for P2 with r = l = h 3 



mesh 


size h 


IHIo 




IMIi 








error 


order 


error 


order 


2 


-l 


6.0137E-01 




2.5073E-00 




2 


-2 


6.3958E-02 


3.2331 


8.5316E-01 


1.5552 


2 


-3 


7.4660E-03 


3.0987 


2.3528E-01 


1.8584 


2 


-4 


9.5200E-04 


2.9713 


6.0522E-02 


1.9588 



We also check the convergence order for the method of characteristics devel- 
oped in Section |3j The corresponding numerical result are provided in Table 
[3] From this table, we can find the convergence order is 1 which is the same 
as in ( 13781) . 

Now we come to check the efficiency of the parallel scheme of Algorithm 15.11 
For this aim, we set the discretization parameters h = 2~ 8 , r = i = 1/512 and 
the linear finite element method is adopted. The consuming time (in seconds) 
are shown in Table [H From Table [U we can find the parallel scheme Algorithm 
15.11 has good expansibility. 

We also check the consuming time in each processor for different scale in 
each processor. For each test, we run 8 time steps (N — 8). Tables [5] and [6] 
show the corresponding consuming time (in seconds) for the average time and 
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Table 3. Error and rate of convergence in the internal coordi- 
nate for the method of characteristics with P2 (h = i) 



mesh 


size h 


ll e llo 




IMIi 








error 


order 


error 


order 


2 


-i 


6.3862E-01 




2.8423E-00 




2 


-3 


3.4562E-01 


0.8858 


1.5382E-00 


0.8858 


2 


-4 


1.7650E-01 


0.9695 


7.8427E-01 


0.9718 


2 


-5 


8.8689E-02 


0.9928 


3.9404E-01 


0.9930 


2 


-6 


4.4398E-02 


0.9983 


1.9726E-01 


0.9980 



Table 4. Strong 


parallel test with PI (h = 2 


8 ), r = 


1/512 


and i = 1/512 








number of processors 


8 16 32 


64 


128 


time (in seconds) 


28103.01 13555.03 6832.26 


3708.71 


1840.43 


rate of speed up 


1.00 2.07 4.11 


7.57 


15.26 



maximum time, respectively, for all the processors. These two tables also show 
that Algorithm 15.11 has good parallel property. 



Table 5. Weak parallel test with PI element (h 
age time in seconds 



aver- 



number in I 


1 


2 


4 


8 


16 


8 


9.30 


15.30 


27.51 


55.28 


116.42 


16 


9.91 


15.44 


28.44 


59.44 


117.24 


32 


9.85 


16.98 


32.02 


60.93 


118.89 


64 


10.01 


17.28 


32.66 


63.88 


121.96 


128 


10.21 


17.98 


33.55 


64.27 


127.63 



Table 6. Weak parallel test with PI element (h — 2 8 ): max- 
imum time in seconds 



number in I 


1 


2 


4 


8 


16 


8 


11.19 


16.10 


27.60 


60.52 


120.26 


16 


11.26 


16.43 


31.54 


61.36 


120.83 


32 


12.73 


18.50 


35.29 


68.18 


131.98 


64 


11.20 


19.63 


36.39 


75.43 


133.55 


128 


12.86 


20.28 


38.01 


73.63 


146.01 



7. Concluding remarks 

In this paper, we are concerned with the parallel numerical method for the 
PBEs with one internal coordinate posed on the domain (0, T] x Q e x fi x with 
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the dimension 1 + 1 + d. The parallel scheme is based on the method of 
characteristics and the finite element discretization. Some numerical results 
are also provided in Section E] to demonstrate the efficiency of the proposed 
method. 

Here, for the simplicity of the description of the numerical method, we as- 
sume the diffusion coefficient e is large enough such that the diffusion is dom- 
inated. For the convection dominated case (c.f. [TJ [TOj [13]), we will combine 
the method of characteristics and the stabilized finite element methods (c.f. 
PQ El ESI ED]) and this is our future work. Furthermore, the parallel method 
should also be applied to the simulation of the industrial crystallization process 
(c.f. |HJE2]) anc ^ °t ner similar models (c.f. [7]). 
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